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In principle, the Luttinger-Ward Green's function formalism allows one to compute simultane- 
ously the total energy and the quasiparticle band structure of a many-body electronic system from 
first principles. We present approximate and exact expressions for the correlation energy within 

, , the GW-RPA approximation that are more amenable to computation and allow for developing 

efficient approximations to the self-energy operator and correlation energy. The exact form is a 
sum over differences between plasmon and interband energies. The approximate forms are based 
on summing over screened interband transitions. We also demonstrate that blind extremization of 
such functionals leads to unphysical results: imposing physical constraints on the allowed solutions 
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(Green's functions) is necessary. Finally, we present some relevant numerical results for atomic 
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systems. 
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1. INTRODUCTION 



Density Functional Theory (DFT) [TJ [2] with the local density (LDA) or generalized gra- 
dient approximations (GGA) [2H3] is the most widely used framework for first principles 
calculations of materials. In practice, it is often found to provide a good description of 
ground-state total energies, atomic geometries, vibrational modes, etc. of a variety of ma- 
terials. A major shortcoming is its inability to predict accurate electronic band structure 
energies [5j. For band insulators with weak correlation the failure is mainly quantitative 
[Bj. However, for the classic case of transition metal oxides, the failures can be qualitative 
and serious such as predicting a metallic instead of an insulating ground state [7] . Ideally, a 
computationally efficient ab initio method with accurate total energies and band structures 
would yield a major advance in predictive power. 

A number of approaches aim to improve electronic band structures. For systems with 
transition metals, one fundamental problem is that the LDA or GGA does not capture the 
proper electronic correlations for the spatially localized d- and /-state derived bands. Two 
current solutions to this deficiency are LDA+U [7\ and dynamical mean-field theory (DMFT) 
[SJ IS]- The LDA+U approach is popular and easy-to-use: one "manually" adds static and 
localized correlation effects within a Hubbard-like model to the DFT energy functional. 
More sophisticated is DMFT where localized but dynamic and high-level correlations are 
included using exact solutions of interacting quantum impurity models. Both approaches 
are physically motivated in that they create frameworks that include to physics deemed 
important for the problem at hand. A shared drawback is their requirement of an unspecified 
localized basis set in order to define key quantities such as the Hubbard parameters or the 
impurity site. This raises questions of transferability, i. e., the dependence of their predictions 
on the chosen orbitals or parameter values. 

A different approach is to use many-body perturbation theory of Green's functions. The 
most successful is Hedin's GW approximation to the electron self-energy ^U\. This approx- 
imation delivers high quality band structures of many band insulators and simple metals 
[BJ [TJ]. The GW approximation includes a great deal of physics including exact exchange, 
localized Coulomb repulsion, dynamic screening, and dispersion forces: e.g., the LDA+U 
is a static approximation to GW [TJ. However, most GW calculations are perturbative in 
that they compute corrections to a mean- field DFT input, and the quality of the final result 
depends on the DFT description. In certain transition metal oxides where the LDA pro- 
vides a decent starting point, GW corrections on this DFT starting point can yield a good 
description of the electronic bands [T2HTTJ. However, in other situations, the inadequacy of 
the DFT description can lead to quantitative errors [121 HSrEO]. 

Clearly, it is advantageous to apply GW beyond the perturbation-off-DFT prescription. 
Such an approach would not assume a localized basis or any set of parameters. Recent 
methods such as the Quasiparticle Self-Consistent GW (QSGW) [IS] or the self-consistent 
COHSEX (scCOHSEX) [20] have successfully moved away from using DFT as the starting 
point. These two methods find the noninteracting initial guess for the band structure ap- 
proximately but self-consistently within GW. Such approaches allow for inclusion of both 
static and dynamic screening effects in addition to some localized (Hubbard U) Coulombic 
physics in a single, general, and parameter-free framework. QS GW and scCOHSEX are 
self-consistent band structure methods, but it would be highly desirable to turn them into 
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total energy methods via the Luttinger-Ward [22] approach that should, in principle, allow 
one to obtain accurate total energies and band structures. 

Separate from self-consistent band structure methods, there has been ongoing work on 
using GW-RPA type correlation functionals for computing total energies. Much of the 
activity was sparked by initial work on the uniform electron gas showing that Luttinger- 
Ward functionals with the GW-RPA correlation energy provided very accurate ground- 
state total energies [23H26] . Other model calculations cast doubt on whether such high 
accuracy was generic (27]. However, actual calculations on atoms, small molecules, and 
simple solids find ground-state energies that improve over the standard DFT functionals, 
especially when short-ranged corrections are added J2EH331 ESH3H] • The vast majority of such 
calculations are post-DFT calculations. Namely, they use Green's functions based on LDA 
or GGA wave functions and eigenvalues to evaluate the correlation energy instead of using 
a self-consistent one-particle description coming from a more elaborate, and presumably 
more accurate, theory like QSGW or scCOHSEX. In addition, the correlation energies are 
evaluated using the standard formula relying on frequency integration of the trace of a 
matrix logarithm (see Section [5]). 

In this work, we report on three main points of progress in this general area for the 
GW-RPA correlation energy functional. The first is an exact rewriting of the correlation 
energy in terms of plasma and interband energies, a result that was found recently using 
very different methods [21J. This exact form is amenable to computations of identical com- 
plexity to present-day linear response time-dependent DFT (TDDFT) or Bethe-Salpeter 
(BSE) calculations [TTI HBTflffi In addition, this new form has much better convergence 
properties when compared to the standard frequency integration method as evinced by our 
atomic calculations. Second, we prove that the GW-RPA correlation energy functional is 
not bounded from below, has a minimum with negative infinite value, and when evaluated 
using non-interacting Green's functions, it has no extrema. This means blind optimization 
of total energy functionals that are based on the GW-RPA correlation energy is highly prob- 
lematic, and physical constraints are required. Third, we rewrite the GW-RPA correlation 
energy approximately as a sum over screened interband transition contributions. This allows 
us to create a ladder of approximations of which the COHSEX pU] is the lowest wrung. In 
addition, this allows us to put schemes such as scCOHSEX on a firm footing by showing 
how they can originate from a variational principle. The ladder provides a series of more 
accurate functionals and associated self-energy operators that may deliver improved band 
structure and Green's functions. We hope that these findings pave the way towards efficient 
and accurate computation of total energies within GW-RPA as well as the creation of effi- 
cient self-consistent schemes for computing total energies and band structures that improve 
over DFT. 

This paper is organized as follows. In Section [2j we describe our notation and provide 
basic definitions. Section [3] provides a brief description and the necessary ingredients of the 
Luttinger-Ward approach. Section|4]describes how the total energy Klein functional we focus 
on in this work actually has no dependence on the choice of non-interacting auxiliary Green's 
function that enters the theory. This greatly simplifies the form of the functional if we restrict 
ourselves to evaluating the total energy on non-interacting Green's functions. We also make 
a connection to the Sham-Schluter equation (UJ HH] appropriate for the variational problem 
associated with the functional. In Section |5j we rewrite the GW-RPA correlation energy 
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in terms of integrals of the standard, time-ordered RPA polarizability. Section [6] provides 
expressions for the derivatives of the total energy functional versus the wave functions and 
eigenenergies or equivalently the static and Hermitian non-local potential determining the 
non-interacting Green's function. In Section [7] we perform an exact rewriting of the GW- 
RPA correlation energy functional in terms of plasma and interband transition contributions. 
Section [8] provides a proof of the unboundedness and lack of extremum of the Klein total 
energy functional when evaluated on non-interacting Green's functions. We discuss what 
this result means in more physical terms, how it relates to other results in the literature, 
and the requirement of constrained optimization that stems from this result. In Section |9j 
we switch gears and derive approximate forms of the GW-RPA correlation energy written in 
terms of contributions from self-interaction-free screened interband transitions. Section [10] 
derives a ladder of approximations that allows us to connect to the COHSEX approximation 
for the self-energy and its associated correlation energy expression. We provide a number 
of approximate correlation energies from which appropriate self-energies can be derived. In 
Section 11, we report numerical results for atoms demonstrating the lack of a lower bound 
to the correlation energy, the superior convergence properties of the exact plasmon form, 
and a tabulation of the quality of the various approximate forms derived in Sections [9] and 



10 We summarize and provide an outlook in Section 12 



2. DEFINITIONS & NOTATION 

In this work, we restrict ourselves to systems with time- in dependent, non-relativistic, 
many-body Hamiltonians so that all response or Green's functions are functions only of 
time differences. Furthermore, we will consider systems with time reversal symmetry so that 
quantities such as Hamiltonians, density matrices, or one-particle states are real-valued. We 
set h = 1 so energies and frequencies are interchangeable. 

Wherever possible, we use matrix notation. For example, the one-particle electron Green's 
function for a time-independent system in the frequency domain, G(x,x',u), is a function 
of three arguments. The x and x' arguments include both spatial coordinates and spin: 
x = (r, a) where r*is a three-vector and a = ±1 labels the two spin projections. In matrix 
notation, we write the matrix G(uj) whose matrix elements are (x\G(u)\x') = G(x,x' ,u). 

The time-ordered, non-interacting, one-particle Green's function Gq(u) is given by 

G (uj) = (col - Ho)- 1 = ^MK. (1) 

— w — e n 

n 

The eigenenergies e n have imaginary parts Ime n = sgn(/i — e„)0 + where + is a positive 
infinitesimal and fi is the Fermi energy. Thus the poles of Go are above the real to axis for 
occupied or valence states, labeled by v so e v < /i, and below the real axis for unoccupied 
or conduction states, labeled by c so e c > \i. The non-interacting one-particle Hamiltonian 
H generates the orthonormal eigenstates |n) and real eigenvalues e n , 

H \n) = e n \n) . 

The wave functions in coordinate space are denoted as ip n (x) = (x\n). In the time domain, 
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we have the time-ordered formula 



Qo(t) = %-^G (u) = - i e(t)Y,\c)(c\e-^ t + te(-t)J2\v)(v\e- ievt - (2) 
The non-interacting density matrix is the standard sum over the occupied states 

/c 
-t 



-oo 27Ti 



We separate a potential U from H via 

Ho = T + f/jon + C/o 

where T is the kinetic operator and U ion the electron-ion interaction potential. The static 
and Hermitian Uo represents approximately the effects of the Coulomb interaction (Hartree, 
exchange, and correlation). For example, in DFT Uq is a local potential that is the sum of 
the Hartree and exchange-correlation potentials. However, in general, we allow for a non- 
local Uq, i.e. Uq{x,x') 7^ for r ^ r'. For a fixed nuclear configuration and thus Ui on , Go is 
determined by Uq and vice versa: 

Go(^)- 1 = ujI-Ho = ujI-T- U ion - U (3) 

The exact, interacting, one-particle Green's function obeys the Dyson equation 

= col — T — U ion -<t) H - ^ xc {u) 

where 4>h is the Hartree potential determined by the electron density n(x) 

and the bare Coulomb operator is 

V(x, x ) 



r _ r n 



The self-energy S xc (a;) is frequency-dependent (dynamic) and non-Hermitian and encodes 
the complex exchange and correlation effects of the many-body system. We can write the 
Dyson equation equally as 

G-\u) = G \u) - [0 + E xc (cj) - Uo] . (4) 

This shows that to obtain the true interacting Green's function, we replace the static, Her- 
mitian Uo by the dynamic, non-Hermitian <pn + £ xc (u;). The interacting electron density 
n(x) and density matrix p(x, x') can be computed from the Green's function via 

/OO J 
— e luj0+ G(x, x', u) = -iQix, x', t = -0+) (5) 
-oo 27r « 
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where the last form is the Green's function in the time domain evaluated for infinitesimal 
negative times. Note that the relations of Eq. ^ among n, p, G and Q also hold for the 
non-interacting n , p , G and Q . 

The frequency dependent dielectric matrix e(u) is related to the irreducible polarizability 
matrix P(oo) via 

e(u) = I - VP(u) . 

Within the RPA approximation, P{oS) is a sum over transitions between valence and con- 
duction states, 

p/ 2(e c - e v )\cv)(cv\ 

c,v ' 

The pair states \cv) are defined in coordinate space via the single-particle wave functions 
ip n (x) through 

(x\cv) = (x\c)(x\v)* = ip c (x)ip v (x)* . 

For a more compact notation, we label each interband transition (c, v) by t with energy 
A t = e c — e v — i0 + and can write 

^) = En ( ( u ) , n t (o;) = ^§, (6) 

where II t (a;) is the polarizability contribution of transition t. We can write P{oS) in a more 
general form to handle possible degeneracies in the transition energies A t . Namely, we break 
the sum over all transitions t into a first sum over the distinct energies A and then a sum 
over all degenerate transitions S with energy A: 

OA 

PM = ^n A ( w ) , n A (c) = ^-^ £ \s)(s\. (7) 

A <5|A 4 =A 

Finally, the time domain RPA polarizability is 

V(x, x', t) = -iG (x, x', t)G (x', x, -t) = [ ^e~ iut P(x, x' , u) . (8) 

3. LUTTINGER-WARD APPROACH 

One of the central points of the overall approach of Luttinger and Ward [22] is that one 
can obtain both the ground-state total energy and interacting G(u) from the extremum of 
an energy functional of G. In this work we concentrate on the specific case of the Klein 
functional [50], a functional of both G and an auxiliary non-interacting Go that is meant to 
be an initial guess for G. It is given by 

F[G, Go] =r ^e iM+ tr\H G (u) + I - Go^G^u) 

+ ln[G (o;)- 1 G(a;)] - U G(u)\ + E H [n) + $ XC [G) . 
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The u integral can be turned into a closed contour integral over the upper complex u plane 
due to the factor e luj0+ . Eh is the Hartree energy stemming from the electron density n(x) 
associated with G(cu): 

Eh\p\ = 2 J d% J dx' n(x)V(x, x')n(x') = - J dx n(x)0#(:r) ■ 

The functional $ m [G] is the exchange-correlation energy functional for this approach and 
is a complicated functional of G. Formally, it is a sum over all diagrams to all orders in 
the Coulomb interaction obtained by closing all skeleton self-energy diagrams with Green's 
functions (with appropriate weight) [22]. Much like in DFT, choosing an approximate form 
for $a; C corresponds to including a certain level of treatment of exchange-correlation effects. 
However, since $ xc is a functional of the more information-rich, non-local and dynamic 
G(x,x',u) as opposed to the simpler density n(x) in DFT, relatively simple forms for & xc 
will be equivalent to rather complex functionals of the density in DFT. On the other hand, 
compared to DFT, the price for more information is the increased complexity of the energy 
functional F and the entire theoretical and computational approach. 

Within Luttinger-Ward theory, one focuses on the extremum of the functional F. At the 
extremum, the value of F is the true ground-state total energy, and the extremizing G is the 
true one-particle Green's function [10, 22j. Therefore this framework provides, in principle, 
both exact total energies and one-particle properties such as quasiparticle wave functions, 
electron densities, and band structures. To find the variation of F versus G, we use a matrix 
differentiation rule based on the properties of the determinant: 

tr{\n{A)} = ln{det(A)} -»■ 8tr{\n(A)} = triA^SA} (9) 
where 6 A is the variation in the matrix A. The variation of F for fixed Go is 

I" ^ e i^o+ tr ( _ Gofay^Gfu) + Giujy^Gioo) - U 6G(u) }+8E H + 8$ xc 

8G(u)\. (10) 



5F 



G 



oo 

00 du 



■e tul0+ tr 

2ni 



Gico)- 1 - Goico)- 1 + 4> u + - U 

dG{U!) 



Setting this to zero for arbitrary 5G yields the Dyson Eq. Q with the self-energy given by 
the functional derivative 

Tj xc (u) = 27T2- 



8G(u) 

Again, the situation is analogous to DFT where the exchange-correlation potential is the 
functional derivative versus the density of the exchange-correlation energy functional. 

Within Luttinger-Ward theory, there are two separate challenges. The first and most ob- 
vious is choosing some approximate form for $ xc . Second, there is the additional challenge 
that, unlike DFT where iV-presentability conditions for the electron density n(x) have been 
known for a long time [SU [52], similar conditions for the Green's function G(x, x', u) are not 
known to the best of our knowledge. In other words, it is not generally known which subset 
of functions G(x, x', u) correspond to physically realizable Greens functions for interacting 
electrons with the standard many-body Hamiltonian. Therefore, one must also decide on 
some scheme for restricting oneself to physically correct forms of G. From a more pragmatic 
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viewpoint, working with an arbitrary function G(x,x',u>) of three variables is computation- 
ally very demanding so that any simplifying assumptions on the form of G are enormously 
helpful in practice. In what follows, we will restrict G{uS) to be of non-interacting form, and 
the next section explains why this is a sensible choice. 



4. LACK OF DEPENDENCE OF F ON G 

Clearly some approximations are required to make progress. We will first replace the 
true interacting G(oS) by something simpler in order to reduce the search space for the 
extremum. Restricting to noninteracting Green's functions that are generated by Hermitian 
Hamiltonians is sensible from a physical viewpoint: we will try to search for the "best" 
noninteracting picture of the electronic system. As an added benefit, we can employ known 
algorithms for Hermitian and orthonormal eigensystems. 

In fact, this restriction greatly simplifies the structure of F. This is because F does not 
in fact depend on the choice of U or equivalently G . Namely, 6F/6U = for fixed G. 
Physically, this is sensible since the final result for G should not depend on the arbitrary 
initial guess Gq. For example, in the Dyson Eq. Q, we remove Uq and replace it by the 
self-energy. To demonstrate this lack of dependence, we use Eq. Q to find the variation of 
F versus G at fixed G: 



5F 



g ./.oo 2-Ki 



<ljJ e lM+ tr{ [8H ]G (u;) + H [5G (uj)] - [5G (a;)- 1 ]G(a;) 



+ G(w) -1 G (o;)[5GoM _1 ]G(u;) - [SU ]G(u) .} 
Using Gq(u) = (ul — Ho)" 1 and the variation of an inverse = — .A -1 [<L4.]A -1 , we have 

5H = 5Uo , dGoiu)- 1 = -5U , 8G (u) = G (u)[8U ]G (u) . (11) 
Plugging these in gives 

8F = r ^e^ 0+ tr{[5U )G {cu) + H G (u)[SU )G (u) + [5U )G(uj) 

- G{u)- l G Q {uS)[5U Q ]G{uS) - [6U ]G(u 
Using the cyclicity of the trace, this simplifies to 

duJ A^ + Jn (, ,m n (, axtt \ - f°° duJ ST tn(n\5U \n) 



SF 



where we have traced over the orthonormal basis of eigenstates |n) and used the diagonal 
nature of Hq and Go i n this basis. The e lul0+ factor has us close the integral over the upper 
oo complex half-plane, 

dio sr-^ e n ( n \8U \n) 



5F 



E 



G J 2iri ^— ' (u) — e n 

n v 



(2 
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The numerators have no u dependence. The standard contour integral for analytic f(z) 



dz f(z) 



d k ~ l f(z) 



2m (z-a) k (k-l)\ dz k ~ x 



(12) 



applied to this case with k = 2 gives a zero derivative and null result. So we have our desired 
result 

SF[G,G ] 



SU 



0. 



G 



(13) 



Hence, for a fixed G, we can evaluate F at any valid G without changing its value. 
Physically, F doesn't depend on the initial choice of non-interacting Hamiltonian. Since we 
are restricting G to be noninteracting, we may as well set G = G . Then F simplifies to 



F[G , G c 



du 
2tH 



e™° tr{ [T + U lon ]G (u) \ + E H [n ] + <S> xc [G oi 



tr|[T + U lon }p } + E H [n } + <$> XC [G C 



(14) 



We have the noninteracting kinetic, electron-ion, Hartree and exchange-correlation energies. 
Operationally, the first three terms are identical to their counterparts in DFT. Due to the 
extremal nature of F about the extremizing G, F[Gq, Go] provides a variational estimate of 
the ground-state energy with the error being smallest for the "best" Gq. Except for $ xc , the 
other energy terms depend only on the density matrix p . Only Q xc depends on the energy 
eigenvalues e n . 

Before we end this section, we make brief comments on the functional derivative of 
F[Gq,Gq] versus Go or equivalently versus the potential Uq. Starting from Eq. (10) with 
G = Gq or Eq. (14) and using the relation of 5Gq to 5Uq in Eq. (11 ), we have at least three 
equivalent ways of writing the variation of F[Gq, Go]: 



SF[Gq, Gq 



-oo 2 ™ 

00 du_ 

2rd 



-oo 

00 du 



2ni 



e lM tr^[<p H + H xc (uj)-Uo}5Go(uj) 
e i « 0+ tr{ [T + U lon + <\> H + 5G (ou)} . 

e^ + tr\Go{uj) [<j> H + E xc (u) - U ] G (lo) 6U 



(15) 
(16) 
(17) 



Eq. (15), written in terms of the variation 5Go, is used below when computing derivatives of 
F versus the ip n {x) and e n that characterize Go- Eq. (17) is written in terms of the variation 
of the potential U that generates Gq. Setting 5F = in Eq. (17) for arbitrary 5U yields 
the matrix equation 



du 



e w0+ Go{u) [<f> H + E xc {u)}G {u;) 



du 



Go{u)UoGo{u) 



(18) 



which is the linearized Sham-Schliiter (LSS) equation for the exchange-correlation potential 
operator V xc = Uq — <pH [471 14*8] . This is the condition determining the static and Hermitian 
V xc that most closely resembles the dynamic and non-Hermitian H xc (u). As Casida has 
noted, this is equivalent to saying that V xc is optimal in a variational sense [53]. The LSS 
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equation is most often discussed in the context of DFT where V xc (x) is a local function 
whereas what we have here is the more general case of a non-local Hermitian V xc (x,x'). 
Solving the LSS is nontrivial because (a) the frequency integral must be evaluated numeri- 
cally, and (b) Go and E xc depend on Uq in a nonlinear manner making for a self-consistent 
problem. We will return to the question of local versus non-local V xc in Section [8] where 
we will see that solving this equation as written generates unphysical results, and we will 
discuss what is known in the literature for local V rr functionals. 



5. IN TERMS OF THE STANDARD RPA POL ARIZ ABILITY 

Beyond using an approximate or constrained form for G, we also require approximations 
for <& xc . In this work, we study <& xc within the GW-RPA approximation, normally defined 
via 

^f[Go] = ^/^tr{ln.»}. (19) 
The dielectric matrix s(u) is modified [lOj and related to a modified polarizability P{ui) via 

e(w) = / - VP(u) . 
The modified P is most easily defined in the time domain 

v{t) = r %-^p(u) 

J-oo ^ 

as 

P(x, x', t) = -iG (x, x', t - + )g (x', x, -t - 0+) . (20) 
Specifically, since p (x,x') = —iQ (x,x',t = — + ), for zero time argument we have 

P(x, x' , 0) = ipo(x, x')po(x', x) . (21) 

This P(t) is modified from the standard P(t) by the infinitesimal shifts in the arguments of 
Go [10]: P{t) and Pit) of Eq. ^ are identical for any finite non-zero t: they only differ in 
an infinitesimal neighborhood about t = 0. This ensures that the Fock exchange energy in 
i s properly recovered [TU]. Hedin has provided formulae for Q xc including corrections 
beyond GW-RPA, so that from a formal point of view, it should be possible to proceed 
beyond the GW-RPA approximation [10J. 



Computing & x ^ from Eq. (19) is cumbersome: we must converge a continuous integral, 
and for each u we need a matrix logarithm. Furthermore, the physical meaning of the 
formula is hard to appreciate, making it difficult to create systematic approximations. One 
of the main aims in this paper is to rewrite ^ x ^ in more tractable forms that permit us to 
understand its physical content. 



However, to proceed, we need to first rewrite Eq. (19) in terms of the usual time-ordered 



dielectric matrix s{u) = I—VP(u) where P(u>) is the standard RPA polarizability of Section 



2l We expand the logarithm in Eq. (19) to all orders using 



x J 



ln(l-*) = -£ , 

3=1 J 
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to arrive at 



The j = 1 term in the series is 



Writing the trace explicitly the x basis and using Eq. (21 ) shows that the j = 1 term is just 
the Fock exchange energy Ex 

j = 1 term : — - J dx J dx' po(x, x')p (x', x)V(x, x') = E x [po] ■ (22) 

Thus, we have 

j=2 J 

This has naturally separated Ex from the correlation energy $ c , 



= $ Tr - E 



x 



The remaining frequency integrals with j > 2 correspond to j th order autocorrelations of 
V(t) evaluated at t = 0, so the infinitesimal time shifts actually do not matter. Specifically, 
the j th term has an u integral 

/oo j /-oo j /-oc />00 

— tr{[VPHp} = / — / ^•••/ dt j e l ^ + - +t ^tr{VP(t 1 )---VP(t J )} 
■oo J —oo J —oo J —oo 

/OO POO 
dh ■ ■ ■ I dtj <J(*i H h i,) tr (VP(ti) • • • VV{tj)} . 
-oo J— oo 

We integrate over the times so that the infinitesimal region about t — where 7^ T 7 ^) 
is of zero measure. Hence, we can safely replace V by V for the j >2 terms to get 

j=2 J J ~°° 

„ r . 1 f°° doo r„ , ,1 I [°° du r , .-I 
= Bx[po] - = 0)} + i jT ^(r{lne( W )' 



In the second line, we added and subtracted the corresponding j = 1 term to the first line 
and summed up the series for the logarithm. In the third line, evaluating V(t = 0) requires 
some care as plugging in t = directly in Eq. (pi) yields an indeterminate result. Rather, 
we take the limit t -> of V{t) in Eq. Q to find 
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or in matrix form 

V(t ^0) = - z J2 \cv)(cv\ =-iJ2 l*X*l • 

c,v t 

The direction of approach of the limit, t — > + versus t — > — + , is immaterial as both yield 
the same answer for real valued Green's functions Go(u) due to time reversal symmetry. We 
then have our desired result: written in terms of the time-ordered s(ou), 

^[G»] = E x [p»]-\Y,(t\V\t) + \ I™ ^-tr{\ns{u)} . (23) 

t J-oo 

This expression is our main workhorse. In Sections [7] and [9j we evaluate this integral exactly 
and approximately to generate alternative forms for $^f c M/ . 



6. DERIVATIVES OF F AND 



In this section, we provide expressions for the derivatives of the exchange-correlation 
functional $^ c w [Go] an d the total energy functional F[G , G ] versus the wave functions ip n 
and eigenenergies e n or equivalently the potential U that determine G . The e n derivatives 
are used in Section [8] to prove the unboundedness of the energy functional. Separately, these 
derivative expressions can prove useful to investigators contemplating variational approaches 
that extremize F[Gq, Go] which require analytical derivatives. 

We begin with variations of the eigenenergies e n . The derivative of Go is 

dG (u) \n)(n\ 



de n (u - e n ) 2 

As per Section |4j the only non-zero contribution to the variation of F[Gq, Go] when changing 
e„ comes from the exchange-correlation functional $^ c w so 

dF[Go,G ] _dQ™[G ] _ r duJ -e^ + trh xc (u) dG ° {Uj) 



de n de n J_ OQ 2ici I <9e 



We turn this into a contour integral over the upper complex u half plane. Using Eq. (12), 
we get contributions from the poles of T, xc (u) that are above the real axis and a possible 
contribution from the pole of dG /de n if n is occupied (if e„ < //). 

To organize the process, we write E xc (o;) in the general form of a sum over poles 

Here, S x is the Fock (bare) exchange operator 

E x (x,x') = - i> v (x)ip v (x')*V(x, x') 
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The are the poles of T, xc (u) that are above the real axis, Im > 0, and <r+ are their 
residues, while and £1 are the residues and poles below the real axis, Im £1 < 0. This 
allows us to write the derivative as 

dF[G ,G ] _ v - (n\a+\n) dZ xc (u) 
A quick manipulation yields 

If the residue matrices a ± are positive definite, then the derivatives for empty states are 
always positive while those of filled states are always negative. For the GW-RPA approxi- 
mation, we can demonstrate this explicitly: as we will discuss in more detail in in Section [7], 
the RPA screened Coulomb interaction W(u) = £ -1 (cu) has the outer product form 

w RpA ( U ) = v+vy 2u F lp)(p } v. 



p p 



A standard exercise in time-dependent perturbation theory shows that the exact many-body 
W also has the same outer product form 

W e,act {u) = V + V Y m-E )\n s0 )(n s0 \ 

where E is the exact ground-state energy, E s are the exact excited-state energies, (x\n s o) = 
(0\n(x)\s) where n(x) is the fermion density operator and |0) and \s) are the many-body 
eigenstates. Thus the exact and RPA W are related via the replacements p, u p , \p) — > 
s, E s — E Q , \n s0 ). 

Whether we have the exact or RPA W, the GW self-energy is 

e~ %w 0+ G(x, x , u — u')W(x, x', us') . 

-oo 2vr 

For G = Gq and the RPA W, the integral yields 

v / / \ v / a , v^/ , T/ | w , T/ | a /V^ ip v (x)ip v (x')* ^ ip c (x)ip c (x'y 

Z xc (x,x,uj) = E x (x,x) + > {x \V\p}(p\V\x) > ■ + > 

/ — ' \ ' cu + u v — t v / — ' cu — w„ — e c 

p \ ti K c K 

Thus the poles and residues of the GW S xc are indexed by an eigenvalue e n and a plasmon 
u p : the labels a or /3 in Eq. (24) correspond a pair index (n,p). The precise correspondence 
for the GW-RPA approximation is 

a+ p (x,x') = (x\V\p)(p\V\x')i) v (x)i) v (x')* , C,p = e v - uj p (26) 

and 

v7, P { x i x ') = <ar|V r |p){p|V'|a/)^ c (ar)V'c(ac')* > Q> = e c + w p • ( 27 ) 
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We arrive at the desired derivative 

dF™'{G a ,G a } v \(nv\V\p)\ 2 „,„ -y \(nc\V\p)\ 2 . . 

^ c,p ' 

It is clear that for unoccupied states the derivative is always positive and for filled states it is 
always negative. Repeating the derivation with the exact W produces the same conclusion 
which is thus inherent to the GW form of the self-energy when G is non-interacting. We 
discuss the implications of this result in Section [8j 

For completeness, we provide the derivatives of F versus the wave functions ip n (x) and 
versus the potential Uq. For the wave function derivatives we have 

5F[G ,G } f°° du 



f 00 ^ + tr{[T + U ion + + E xc (u)]^fl] 

J-oo 2 ™ 1 Ol/j n {x)*J 



The functional derivative of Go is 

SGq(x', x", Uj) 1pn(x ; )5(x" — x) 



5ip n (x)* oj-e n - zO+sgn(/i - e n ) ' 
We insert this into the integral and perform the contour integral. For the general case, 

5F J?°;t° ] = 8(» - tn)(x\T + U ion + ^ H + E xc (e n )\n) + £ (29) 
= 9(fi - e n )(x\T + U ion + 4>h + E x \n) 

a Set £n a ^ct € n 

and expressions for the GW-KPA approximation are found by performing the substitutions 
of Eqs. ( 26|27 ). For the derivative versus Uq, it is convenient to express the potential in the 
one-particle eigenbasis so U j k = (j\U \k). The variation 5U is then 

SU = Y / \j)6U 0jtk (k\. 

3,k 



Inserting this into Eq. (17) and performing the contour integral and some algebraic manip- 
ulations yields 

v\(j) H + E xc (e v ) -U \c) ^ (vWM (o \ 



dF[G , 


G ] 


dU 0v 


v> 


dF[G , 




dU 0cv 


dF[G , 


Go] 




,c 


dF[G , 


Go] 



(c\(j)H + Z xc (e v ) - U \v) (cWHv) 



V K - 1 al 1 > (33) 

y . (34 ) 
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where v, v' label valence states and c, d label conduction states. As a check, the diagonal 
cases v — v' or c = d yield the same result as the e n derivative as necessary by first order 
perturbation theory. Setting all the derivatives to zero is equivalent to solving the LSS 
Eq. (17) with no constraints on U . This turns out to yield unphysical results as explained 
in Section HI 



7. REWRITTEN EXACTLY WITH PLASMON AND INTERBAND ENER- 

GIES 

In this section, we rewrite exactly in terms of a sum over plasmon and interband 

transition energies. This exact expression has been reported by other workers using different 
methods [21]. (We found this result contemporaneously and independently.) Our derivation 
is based on contour methods which we present for completeness before moving on to the 
implications of this exact expression. 



We start with the result of Eq. (23). The time-ordered RPA s{u) is given by 

2A t \t)(t\ 



e(u) = 1- VP(u) =I-VJ2 



u 2 - A? 



which has poles at the interband energies u = A t . The poles of the inverse dielectric 
function e{oo)~ l occur at the plasma frequencies which are the natural modes of collective 
charge oscillation. The matrix ^(u;) -1 is given by 



e(u )-. = , + V X ( U ) = / + Ki:£!M (35) 



(The uj p have infinitesimal negative imaginary parts.) The plasma frequencies u p and mode 
functions \p) obey the RPA/Casida Hermitian eigenvalue equations 



&tCt, P + Y, 2 V^~'(t\V\t')C t , tP = co 2 p C ttP . (36) 
t 

The column vectors Ct tP are orthonormal and the \p) are given by 

|p) = £i*)J-a P . 
t V u p 

We see that the u p are the eigenvalues of the square root of the Q 2 matrix defined as 

fi? )t , = A t V + 2y/A t AAt\V\t') . (37) 

Before going forward, we make a clarification concerning our use of the terms "plasmon" 
or "plasma" modes or frequencies. In this work, these terms refer to solutions of the 



Casida/RPA Eq. (36) so that the u p are real frequencies. Of course, a plasmon excita- 



tion in an actual material has finite life time due to damping processes. Thus, we are 
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dealing with the real part of the plasma frequencies, and this approximation is analogous 
to our assumption of a non-interacting Green's function where the quasiparticles have real 
eigenenergies and infinite lifetimes. 

Since the poles of s{u)~ l are at u p , its inverse e(u) has a zero at the u p , i.e. e{u p ) has a 
zero eigenvalue and is not invertible. Therefore, for u = A t one of the eigenvalues of e(u) 
diverges while for u = u p one of its eigenvalues is zero. The number of transitions t and 



plasmons p are equal because Eq. (36) is a square matrix Hermitian eigenproblem. 



Letting the eigenvalues of e(u) be \ m (u), the trace over the logarithm of e(u) is 
*™ [Go] = E x [p ] - \ 5>|V|f> + \ ^ $> A m (c) . 

We close the integral on the complex u lower half plane to yield a closed contour integral. 
The function In z is analytic everywhere except on the branch cut starting at z = and 
extending to z = — 00 along the negative real axis. This means that the integrand is analytic 
everywhere except between the zeros and poles of X m (u), i.e. between some interband energy 
A t where 

, , a \ constant . , n . 

\ m (u ->■ A t ) = — + 0{{u - A t )°) 

u — l\t 

and some plasma frequency u p where 

\ m (uj — > Up) = constant • (co — u p ) + 0((oj — co p ) 2 ) ■ 



Due to the positive definite nature of the Coulomb interaction matrix (i|V|i') in Eq. (36), 
the smallest A t is smaller than the smallest u p while the largest u p is larger than the largest 
A t . Hence, we have a set of overlapping branch cuts beginning at the A t and terminating 
at the u p . The contour integral collapses to an integral around the finite segments of length 
Up — A t . We perform the contour integral by remembering that In z changes its imaginary 
part by 2m when going from above to below the branch cut. Each integral about a segment 
yields precisely u p — A t . This leads to our first central result which is an exact rewriting of 
$^ c w/ in terms of plasma and interband energies: 



GW 



xc 



[Go] = E x [p ) + + <*i y l*>) • ( 38 ) 



Since we sum over all t and p, the precise pairings of A t and u p for each branch cut are 
irrelevant for the final result. 

This exact form is much in the spirit of early work on the electron gas showing that 
the RPA correlation energy represents replacing interband oscillators by plasma modes |JT] • 
We refer the interested reader to Ref. [21] for a more detailed discussion of this and related 
points. From a pragmatic viewpoint, this expression is amenable to direct computation since 
the RPA/Casida or related equations are solved regularly when performing TDDFT or BSE 



calculations [HI H9]. Direct diagonalization of the RPA/Casida Eq. (36) makes for 0(N 6 ) 
scaling for a system of iV atoms. However, since only the trace over all modes is needed, 
there are possibilities for improved scaling. For example, the sum over plasma frequencies 
is also the trace of the square root of the Q 2 matrix, 



E 



up = tr { (fi 2 ) 1/2 | , (39) 
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so efficient matrix square roots algorithms that avoid diagonalization are computationally 
advantageous |2T] . 



In terms of convergence, Eq. (38) is superior to the standard integral form of Eqs. (19) 
or (23). There is only a single convergence parameter which is the size of the basis set of 
transitions {\t)} in the calculation. For the integral forms, one must converge additionally 
the continuous oj integral by dealing with issues of grid spacings, truncation at large u, etc. 
In Section 1 1 , we will highlight the superior convergence properties of Eq. ( 38 ) numerically 
for atoms. 



8. UNBOUNDEDNESS OF $£f c w AND LACK OF EXTREMUM OF F 

In this section, we prove that [Go] is in fact unbounded from below: its minimum 
value is negative infinity. Furthermore, the total energy functional F can not have an ex- 
tremum when optimized freely over non-interacting Green's functions Go that are generated 
by non-local potentials Uq. We then discuss the meaning and consequences of these facts. 



Beginning with the total energy functional of Eq. (14), we see that the only term depend- 



ing on the e n is the exchange- correlation part $ xc [Go]- For the GIF-RPA approximation, 



we have the exact expression of Eq. (38) for $^ C H/ . This functional is unbounded from 
below if varied freely over e n . Consider sending all the eigenvalues towards the Fermi en- 
ergy, e n — > p. This makes all transition energies tend to zero from above, A t — > + . The 



RPA/Casida Eq. (36) immediately shows that the plasmon energies must also tend to zero 
Up — » + as well. Thus $XT will turn into 



9™ [G ] -> E x [po] - \ 5>|V|f> = E x [ Po ] - \ EM^M • 

t V c 

Writing this out explicitly in coordinate space 

^[Go]^E x [ Po }-^J2J2 I dx I dx ' ^c{xT^{x)V{x,x')^ c {x')^ v {x'T . 

V c 

Using completeness 

^2ip c (x)ilj c (x')* = 6(x - x') - y]ip v >(x)ip V !(x')* 

c v' 

where v' labels occupied states, we find 

*2HG ] E x [p ] - \Y,j dx Mx)V(x,x)Mx) + \ Y,^v'\V\vv') . 

v v,v' 

Since V(x,x') = 5 aia i/\r — r'\ diverges to positive infinity for x = x', we conclude that 
diverges to negative infinity when all eigenvalues approach the Fermi energy. This 
particular choice of variation of the e n is perfectly permissible if one allows for arbitrary 
variations of the non-local potential Uq. We give numerical evidence of this behavior for 
atoms in Section [Til 
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Therefore, global minimization of F[Gq, Go] over all allowed Go — more precisely over 
all non-local potentials Uq that generate the Go — is guaranteed to give unphysical results. 
However, one may still hope to find a local minimum or an extremum: after all, the Luttinger- 
Ward approach only guarantees that the true G is an extremum of the energy functional. 
However, this is also not possible because dF[G , G ]/de n is never zero. Referring back to 



Eq. ( 28 ) in Section |6j we see that the derivative versus e n is always positive for unoccupied 
states and always negative for occupied states. Hence, a search for an extremum also drives 
us to the unphysical global minimum. No other extremum exists because the derivatives are 
never zero. This proof confirms the suggestion [51] that direct minimization of the GW-KPA 
total energy will give meaningless results. 

We now discuss the meaning and relation of this result to what is known in the literature. 
A number of studies have used local (in coordinate x) exchange-correlation potentials V xc = 
Uo — <Ph when generating Go for atoms and molecules. The choices have included LDA, GGA, 
or optimized effective potentials In addition, a more recent work [37] has solved the 

linearized Sham-Schulter (LSS) equation for spherical atoms to find the optimal local V xc . 
No sign of any instability has been found in any of these works and the numerical evidence 
shows that the LSS-optimal local V xc locates an energy minimum [37]. In addition, there 
are strong theoretical arguments for why the total energy functional will have a minimum 
when varied over local potentials [35J. 

Unfortunately, the optimal local potential is not necessarily the best potential in terms 
of physical predictions. In practice, using a non-local potential to generate G can greatly 
improve results: e.g., for atoms and diatomic molecules, the Hartree-Fock non-local potential 
yields total energies that are essentially indistinguishable from those given by using the 
exact, self-consistent Green's function G [55]. Therefore, including some nonlocality in 
the potential Uq is important physically: for atoms and small molecules, the Hartree-Fock 
description is significantly closer to the true Green's function than any description based 
on a local potential [35]. This result is sensible since the true Green's function G obeying 
Dyson's Eq. Q is generated by a non-local self-energy operator T, xc (u}), and we would expect 
a non-local Uq to be a better approximation than a local operator. 

Alas, our proof above shows that choosing the proper form of nonlocality is not straight- 
forward. The naive idea of extremizing the total energy functional over all non-local Uq 
to find the "best" Go leads to unphysical results. For example, the simple and appealing 
idea of holding the wave functions ij) n (x) fixed and varying the eigenenergies e n to find an 
improved band structure via optimizing the total energy leads to the unphysical minimum. 
Of course, we know that there is an exact G that extremizes the total energy functional and 
solves the Dyson equation self-consistently: there are examples of such calculations in the 
literature for model systems such as atoms, small molecules, or the electron gas [23H26], 135] . 
The problem is that the energy functional has no extremum in the subspace of Green's 
functions that correspond to non-interacting Green's functions. To put this in pictures, the 
simplest likely scenario is illustrated schematically in Figure [T] where one assumes that the 
total energy functional has a single extremum corresponding to the true Green's function. 
Obviously, that extremum must occur for a dynamic and non-local self-energy H xc (x,x',u). 

The situation here is quite different from DFT where one can represent the electron 
density in terms of occupied non-interacting wave functions and where an unconstrained 
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Z(q))-U. 




FIG. 1: Schematic figure showing the simplest likely scenario for the Klein total energy functional. 
What is shown are level curves of the total energy functional F[G] as a function of the self- 
energy Ti xc {uj) that determines the Green's function G via the Dyson Eq. Q (not necessarily self- 
consistently). The horizontal axis represents self-energies that are static and Hermitian operators 
U$ that will generate non-interacting Green's functions Go via Eq. ([3]). The vertical axis represents 
the deviation of the self-energy from a static and Hermitian operator. The black circle represents 
the extremum of the total energy functional corresponding to the true self-consistent self-energy 
which must occur for a dynamic and non-Hermitian Ti xc {uj) since F has no extremum along the 
horizontal axis. 

minimization over allowed densities, or equivalently over any set of occupied orthonormal 
non-interacting wave functions, leads to a single minimum with the correct ground-state 
energy [TJ [2] . For the Klein functional, the analogous approach fails completely because an 
unconstrained optimization over non-interacting Green's functions drives the system to an 
unphysical minimum with negative infinite energy. The analogous situation in DFT would 
be the (fictitious) situation where the correct ground-state density locates a minimum but 
the total energy functional has no minimum or lower bound when evaluated on a subset of 
allowed densities. 

In our mind, there are two basic ways to overcome this hurdle. One could decide to 
solve the Dyson equation for the true self-consistent G and avoid the instabilities. However, 
this requires storing the entire Green's function G(x,x',u) — a full matrix as a function of 
continuous frequency u — which is highly prohibitive in terms of storage and computation 
and unlikely to lead to a method that will deal with realistic systems with many atoms in 
the near future. The other approach is to stick with the simple and appealing picture of a 
noninteracting Go generated by a static and non-local Uq but to constrain the extremization 
so as to avoid unphysical behavior. In other words, one constrains the allowed forms of Uo 
in some way. In this light, the use of strictly local V xc , the QSGW, and the scCOHSEX 
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approaches can be viewed as schemes where one manually imposes physical constraints on 
the potential Uq in order to avoid pathologies. As noted above, when constraining Uq to be 
a local potential, the available evidence strongly argues that the Klein functional will have 
a minimum [33l |35j [37]. However, the open question is how to improve beyond these choices 
and to find better constraints that still avoid pathologies. This will involve coming up with 
some type of presently unknown metric to rank the various types of constraints. 



9. APPROXIMATE ^S 47 BASED ON THE SCREENED INTERACTION 



Although the rewriting represented by Eq. (38) is exact, there are a number of reasons to 



look for other expressions. Ideally, it would be nice to rewrite m terms of the central 
quantity in GW which is the screened Coulomb interaction W . One reason is to remove the 
dependence on the plasmon description and instead to deal only with interband transitions 
and screening. Another reason is that dielectric functions and screening are relatively well 
understood and studied objects for which a variety of approximations and computational 



approaches exist. As we show in Section 10 when deriving the COHSEX and related energy 
expressions, working with the screened interaction allows one to naturally build in different 
types of physical insights. Before wading into the derivation, we highlight the endpoint: 



they key result is Eq. (44) 



We begin with Eq. (23) and re-expand the logarithm, 

*2" [Go] = E x [p ] -If; 1 T^. tr{[VP(u)}>} . 



J=2 



We concentrate on the correlation part $ c = $ xc — Ex and write each P(u) as a sum over 
transition energies as per Eq. ([7]): 

*?"' = - \ E t f t * { E E nw«) •••E ^ m } ■ 

Z j=2 3 J ~°° Z7Ct A x A 2 A,- 

We define Ij as the j th term in this series, 

i s = -^ r^fr{E vn ^( w )E^( w )---E^>)}- ( 4 °) 

•* Ai A 2 Aj 

We close the contour integral over the lower complex u half-plane (the final result is un- 
changed if we choose the upper half-plane). Since each factor Ha(u) is given by 

OA 

n ^ = (.-A) ( , + A) Eim 

the integrand for Ij is analytic everywhere in the lower half-plane except when u = A for 
some transition energy A. Since II a {uj) diverges as (u> — A) -1 close to such a point, the 
integrand for Ij has divergences of type (u — A)~ k where 1 < k < j and k labels how many 
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of the {Ai, . . . , Aj} happen to be coincident. This leads us to write Ij as a sum over sets of 
contributions 7* labeled by the divergence index k, 



j 

k 



k=l 

where 

A qo,—,qk 



1" = - r 

3 2jJ 



/t^E E '^{[^A(o;)] 9o ynA(a;)[FPA(a;)] 9l ---ynA(a;)[FPA(a;)]^}. (41) 



The contour integral is along the real axis and closed over the lower complex u half-plane, 
P a (oj) is the polarization with transition energy A missing, 

Pa(oo) = p(u) - n A (w) = Yl nA '( w ) ' 

and {g , ■ ■ ■ , g^} are any nonnegative integers restricted to sum to j — k, 

k 
1=0 



The prime over the second sum denotes this restriction on {g , ■ ■ ■ , g^}. Eq. (41) is to be 
understood as follows: to get a contribution to Ij, we must have k of the j transition 
energies in Eq. (40) have the same value which we call A while all the remaining j — k 
transition energies must be different from A; most generally, the first go factors in Eq. (40) 
have transitions differing from A, then there is a transition at A, then qi follow which differ, 
followed by another at A, etc.; summing over all A and all possible {qi} includes all the 
possibilities. 

Using the cyclicity of the trace, we combine the first go an d last g& terms, 



J * = ~Tj j E E 'tr{vU A (co)[VP A (co)rVU A (u) ■ ■ ■ VU A (co)[VP A (co)r + ^ 



A go,— m 

and then replace the IIa by explicit sums over transitions to get 

t* = L v (2A)fc v ' v 

3 2j J 2m V - A 2 ) fc ^ ^ 

J J A 1 q ,...,q k Si,-,S k 

tr{y|5 1 )(5 1 i[yP A (o;)]^F|5 2 }(5 2 |---y|4}(4|[^PA(a;)] 90+9& 

= -h <f ^ E J 2 _ A L E ' E (*ii^A( W rw...<(5 Jfc |[7P A ( W )r^y| < y 1 > 

J J A ' q ,...,q k <?!,••• A 

where the transition • • • , ^} sum only over those with energy A. 

In the above expression for ij% there is no separate dependence on g or q^ but only on 
their sum r = g + g^. There are r + 1 possibilities for g and g^ at fixed r. Summing over 
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them generates a multiplicative factor of r + 1. We rename r back to q k and have 

r* = _l / ^ y (2A) fc v ' v 
J 2j / 2ni ^ (cu 2 - A 2 ) k ^ ^ 

J A ' qi,...,q k 

(5l|[^AH] 9i y|5 2 )(5 2 |[yPAH] 9i y|5 3 > • • • + l)(5fc|[^PA(w)] 9fc F|(Jl) • 

We sum over one fewer {<^}, namely from q 1 to g^, where the prime indicates that q± + - — h 
q k — j — k as before. 

The qi are not treated symmetrically in the above expression: qj, is singled out by the 
extra factor q k + 1 because we chose to eliminate go- However, the final result for I 1 - has 
the same value if instead we single out another qi to have the qi + 1 factor because we 
can rearrange the k multiplied factors and permute the q\ via relabeling. Therefore, we 
symmetrize by summing over all k choices of qi (being singled out) and dividing by k. This 
amounts to averaging the qi + 1 factor: 

£E(« + i) = ^-* + *) = { 

i=i 

which happily cancels the 1/j factor. So we now have the symmetric expression 

/ j = 4/^^ fc-f^-A^ £ (5i|[^A(u;)]^|5 2 )...^|[yp A H]-y|5 1 ), 

A V ; <?i,.-,<?fc S 1: -,S k 

We are now ready to sum over j to eliminate the restriction over the {qi}, i.e. remove the 
prime. To do this, we reorder the j and k sums in the correlation energy 

oo j oo oo 

j=2 k=l k=l j=k 

where we added and subtracted I\ = \J2 t {t\V\t). The inner sum over j removes the 
constraint over the qi so 

£>-^£m&f £ £ wm,)-wp.rv\h). 

j=k A Qi,-,QkSi,-,S k 

Each factor of \VP&)} qi V can be summed separately to yield an identical result, and in each 
case we are summing the geometric series for (1 — x)^ 1 : 

oo 

W A (u) = Y^[VPA(ou)] qi V = {I- VP±{u))- l V . (42) 
m=o 

We have defined the modified screened interaction W&{u) for which the interband transition 
energy A is missing from the screening action. Thus the sum over j has yielded 

K-^Ef^ £ www-www. 

j=k A v > S u -,S k 
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We define a square matrix W a {oj) that contains the matrix elements of W A (u) among all 
the degenerate transitions {Si, • • • , 5k} of energy A: 



W A (u) M = (5 p \W A (oj)\5 q ), 
which allows us to compactify the above expression as 

du (2A) fe 



E 7 

j=k 



2-ni ^ k ■ {uj 2 - A 2 ) k 



(43) 



The result of Eq. (43) is exact and summing it over k will recover the same answer as 



the exact result in Eq. (38) albeit expressed in a completely different manner. We can, in 
principle, perform the contour integral. On the lower complex uj half-plane, the integrand 
has a k th order pole at uj = A as well as a large number of other poles coming from W a {uj) 
which physically correspond to the screening modes or plasma frequencies. We would have 
to sum over all residues to obtain the integral exactly. However, the resulting expression is 
unwieldy and can not be simplified in any meaningful manner known to us. 

Therefore, we make a physical approximation to make the integral tractable. The basic 
approximation is in the spirit of the COHSEX approximation: we assume that the physically 
important plasma frequencies are at much higher energies than the dominant interband 
transition energies. This is generally reasonable for solids and extended systems where 
plasmons are strongly collective modes with high frequencies due to the long range of the 
Coulomb interaction. In other words, we assume that the screening dynamics is much faster 
than the interband dynamics. We return to this point at the end of Section [10] where we 
discuss what type of physics is and is not included in this type of approximation. Separately, 
we provide some evidence of the relatively good quality of this type of approximation for 



atomic systems in Section [TT 

Mathematically, this approximation means that we ignore the contributions of the 
residues coming from W A (oS) itself and instead include only the residue of the low-energy 
pole at A. If uj is the energy of a typical pole of W A [uj\ then the neglected terms are 
proportional to powers of the dimensionless ratio A/a). Therefore, if we write the integral as 
a power series in A/a), our approximation amounts to keep only the leading terms of order 
(A/a;) . Equivalently, we pretend that W a (uj) is smooth and analytic so that the only poles 
in the integral of Eq. (43) come from the denominator at o> = A. 

Since W a (uj) is built from the polarization P a (uj), both are missing transitions at energy 
A and are well behaved at and about o> = A. We use Eq. (12) to find 



E 7 

j=k 



k ^ 1 v (2A) k d k ^ 
3 ' ~ 2 4^ k\ duj^ 1 



A) k 



Again, due to the assumption of the smoothness of W A (cu) at low frequencies, its derivatives 
at uj = A are also assumed negligible compared to the derivatives of (a> + A)~ k : this amounts 
to discarding terms with positive powers of A/a). So we arrive at 



E 7 . 

j=k 



J ' ~ 2 ^ k\ 



W A (A) 



d 



k-l 



doj k 



1 



{uj + Af 
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Taking the k — 1 derivatives and evaluating at u> = A yields 



l^(2A) fc ^ f ( — /AN \fc -| (-l) fc " 1 (2A:-2)! 

(2A) 2fc - 1 (A;- 1)! 



2^ 2A J (*!).. (1-2*)/ 

Interestingly, we recognize this as a term in the Taylor series for y/1 + x, 

Vi + X 2^ J („!)2(i_2n)- 

n=0 V 1 V ' 

We sum our approximate expression over all k to get 



oo oo 



k=l j=k 




\Y1 tr { V /A2 + 2AW/ a(A) - /A I . 



2 

A 

where a matrix square root are understood. Putting this together with Jj, we have our main 
result for the approximate rewriting of the correlation energy: 

®c W ^\j2 tr { [ IA2 + 2AWa(A)] 1/2 } - \ (A* + <t| V|f>) . (44) 

A t 

We have achieved our objective of summing the entire series for the GW-RPA correlation 
energy approximately and writing it as the sum of expectations of a screened interaction 
over interband transitions. The modified screened interaction is built from a modified 
polarizability which includes polarization contributions from of all interband electron- 
hole fluctuations except for those at energy A. This means the expression is self-interaction 
corrected in that excited electron-hole pairs at energy A do not screen themselves but are 
only screened by the other transitions. 

Another perspective on this self-interaction correction is provided by identifying con- 



tributions to Eq. (44) that are exact and do not depend on the high-frequency screening 
approximation. These contributions are the terms with j = k (ij terms) because these terms 
have no VP& factors since q — • • • — q% = is forced. Therefore, only the bare Coulomb 
interaction is relevant for them and they represent a weak screening limit. If we include only 
these terms, the high-frequency approximation is unnecessary since all W&(u) are replaced 



by V, and Eq. (44) becomes 



®c W -+ \ J2tr{ [IA 2 + 2AF] 1/2 }-\Y.( At + <W>) • 

A t 

The matrices under the square roots are precisely the block-diagonals of the RPA/Casida 



Q matrix of Eq. (37) with degenerate transition energies. Namely, this formula is the 
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expression we would obtain if we had taken the exact expression of Eq. (38) and computed 



by only solving the block diagonal parts of the RPA/Casida Eq. (36). This provides 



a different viewpoint on why the screened interaction Wa appearing in Eq. (44) must be 



based on a polarizability that does not include any transitions at energy A: interactions of 
transitions with energy A among themselves are already included exactly by this formula 
with no screening (i.e. when Wa = V). So the screened interaction in Eq. (44) must 



be responsible for capturing the couplings among transitions of differing energies which is 
precisely why transitions at A do not contribute to Wa- 

At this point, a great deal of simplification is achieved if the only degeneracies present are 
normal ones (i.e. those due to a symmetry such as spin degeneracy for an spin unpolarized 
system or molecular or crystalline symmetry in real space). Generically, we expect this to 
be the case as any weak perturbation will remove an accidental degeneracy. For the case of 
normal degeneracy, the degenerate subspace spanned by the transitions {S p } with energy A 
transforms as an irreducible representation of the symmetry group. On the other hand, the 
screening matrix Wa(w) must transform as the identical representation (i.e., as a scalar). 
Thus the only non-zero entries in the matrix W&(uj) will be the diagonal entries, and the 
diagonals will all be equal by symmetry. The matrix square root is trivial in this diagonal 



basis and Eq. (44) becomes 



GW 



In 



A? 



2A t (t\W At (A t )\t)-A t -(t\V\t). 



(45) 



This relation can be viewed as a scalar version of the more general matrix Eq. (44). In the 



unusual case of accidental degeneracies, the matrix Wa(u) will have off-diagonal entries and 
it is necessary to evaluate the matrix square root of Eq. (|44]) . Since these degeneracies are 
rare, the matrices in question will be small and using any method (such as diagonalization) 
will not impact the computational load. 

We note that these approximate forms of the GW-RPA correlation energy also suffer 
from the unboundedness problem discussed in Section [8] that plagues the exact GW-RPA 
correlation energy of Eq. (38). For example, consider scaling all interband transitions by a 



factor A, so At — > A At, and then sending A to zero. As long as matrix elements of Wa remain 
finite, the correlation energies of Eqs. (44|45) both diverge to negative infinity for the same 
reason that Eq. (38) diverges. (In fact, when all interband transition energies are scaled to 
zero, the system will show very effective metallic screening and the matrix elements of Wa 
will actually go to zero.) We present numerical evidence of this divergence in Section 11 



We end this section with some comments about the actual usage of Eqs. d44[45[ ). Com- 
putationally, using these expressions requires re-computation of the screening for each tran- 
sition energy. In a naive implementation, the full matrix is recomputed, and since each 
computation of P scales as 0(N 4 ), the overall computational load for evaluating Eq. (44) 



is 0(N 6 ) and thus no better than the exact diagonalization approach of Eqs. (38) and (36) 



Furthermore, the sums in Eqs. p4p5] ) are very difficult to converge numerically because 
transition energies in the continuum are very closely spaced and thus Wa(u) will have con- 
tributions from many poles near u — A. Converging such a sum requires a very dense 
representation of the continuum together with some regularization procedure to control the 



very large contributions from nearby poles. In Section 11 we discuss this problem for atomic 
systems. Thus, the utility of Eqs. (44|45) is not for numerical computation per se but rather 
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for deriving new approximations, as shown in the next section. 



10. COHSEX-TYPE CORRELATION FUNCTION ALS 



In addition to being closed form expressions for the correlation energy, Eqs. (44) and (45) 
are based on the screened interaction. As explained above, literal implementation of these 
formulae is computationally expensive and numerically difficult to converge. Their utility, 
rather, is in analytical work where one can contemplate a variety of approximations to the 
screening that incorporate different physical effects. As an example, we present the simplest 
approximation that leads to COHSEX and its associated <J> C . During the process, we will 
derive a ladder of related approximations. Looking ahead, the relative quality of a number 



of the approximations and simplifications are tested numerically on atoms in Section 11 



We base our derivation on the scalar expression of Eq. (45) for simplicity of presentation. 



The full matrix expression of Eq. (44) can also be used but produces more complex-looking 



results with the same physical content. In what follows, we use Wt(u) as a shorthand for 
W At (u). 

The first approximation is to work within the basic idea of COHSEX: replace the dynamic 
screening matrix by a static one at u — 0. This leads to a first approximate form, 



*f = i^A? + 2A t (t|W t (0)|t> — A t — (t\V\t) 



(46) 



There is a simple relation between matrix elements of the usual screened interaction W(u) = 
£ _1 (w)V A and the W t (oj) appearing here that we exploit to compute the static screening. The 



following matrix identity 



xuv)) 1 u 



u 



-1 



■ft 



/(l - u ] A 1 ux) 



for vector u, matrix A, and scalar x, allows to make the connection. With A = V 1 — P& t , 
u = \t) and x = 2A t /(u 2 



A 2 ), we find 



(t\W t (u)\t) 
For static screening, we have 



(t\W(u)\t) 



1 + (t\W(u)\t) ■ 2A t (u 2 - A 2 )- 1 ' 



(t\W t (0)\t) 



(t\W(0)\t) 



2{t\W(0)\t)/A t 



(47) 



(48) 



Therefore, in terms of computational complexity, evaluation of Eq. (46) requires one calcu- 



lation of the matrix VT(0), computation of its diagonal elements in the \t) basis, and use of 
Eq. (48). For a system of N atoms, this entire project scales as 0(N 4 ). For comparison, 



the original expression of Eq. (44) requires evaluation of W t (A t ) at each interband energy 
A t and scales as 0(N 6 ) in a naive implementation. 



The next approximation expands the square root in Eq. (46) in powers of (t\W t (0)\t) / A t 
and keeps the lowest order term, 



Qstat = 1 J2(t\W t (0) - V\t) + 0((W) 2 /A) . 



(49) 
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This approximation requires that |(£|Wt(0)|i)/At| <C 1. Section [TT] provides numerical results 
for atoms where we show that this expansion of the square root is actually quite accurate. 
Since screening is rather weak in the localized atomic limit, we expect this approximation 
to work even better in extended solids which provide much more effective screening of the 
Coulomb interaction. 

The third approximation is to replace the matrix element (t|Wt(0)|t) by (t|W(0)|t). This 
also creates an error of 0((W) 2 / A) and should be applicable in the same cases as the 
previous approximation. The result is 



= _ ^2(t\W{0) - V\t) + o((W) 2 /A 



(50) 



Completeness \ip c )(ip c \ = I — Y^, v \4 > v){^v\ permits us to rewrite & s c tat using only valence 
states. Putting Ex back and dropping terms of order 0((W) 2 / A), 



stat 



dx p (x, x) W pol (x, x, 0) 



dx / dx'\p (x,x')\ 2 W{x,x',0) (51) 



where W po1 = W — V is the induced (polarization) part of W. The first term is the Coulomb 
hole energy (COH) and the second term the screened exchange energy (SEX). Differentiating 
versus the non-interacting density matrix po yields the static self-energy E^ c at = S^^/Spo, 



V£*M = -W* el (x,x,0)6(x 



x') - p (x, x')W(x, x', 0) + O(8W/5p ) 



(52) 



This is the COHSEX self-energy if we drop the derivative term 5W/5pQ. Ignoring the 
derivatives underlies the successful BSE approach for optical excitations jl6] and excited 
states forces [55]: the successes of the BSE approach suggests that these terms are not 
significant in practice. In brief, the COHSEX self-energy of Eq. (52) is approximately 
associated with and derived from the exchange-correlation energy of Eq. (51). The key 



approximations underlying COHSEX are the static screening limit and the assumption that 
the matrix element (i|W(0)|i) for a transition t is much smaller than its energy A f . Given 
W, Eqs. (51) and (52) depend only on the valence states through p . Having the energy 
functional available in Eq. (51), we also know what missing terms we must add to the 



COHSEX self-energy of Eq. (52) in order to make it a variational scheme. 



To improve COHSEX, we step backwards through the approximations. The first stop 



is to use Eq. (46): this is still a static approximation but doesn't assume that (W)/A is 
small. To go beyond Eq. (46), we need dynamic screening. In order to avoid the 0(N 6 ) 
scaling of the direct implementation of Eq. (44), we can contemplate the following scheme: 
we relate W t (A t ) by to W(A t ) via Eq. (47) and then approximate W by using sum rules 
as per plasmon-pole models to approximate the frequency dependence of W{u) 0, 1561 - 
[58] . Alternatively, we can use model dielectric screening functions [59H6T] to construct 



approximate W(u) and then use Eq. (47) to find W t . Regardless of the specifics, having a 
total energy expression allows us to find the associated self-energy through differentiation. 

We conclude this section with some observations on what is and is not included when 
using the approximate expressions derived in this and the previous section. We expect these 
approximations to certainly include static and localized Coulombic effects automatically. 
This includes screening effects of the medium as well as localized Coulombic physics: even 
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the simplest COHSEX self-energy of Eq. (52) and COHSEX correlation energy of Eq. (51) 
contain the non-local density matrix p that projects the action of the screened Coulomb 
interaction W onto the occupied states. If the occupied states are derived from localized 
states such as 3d or 4/ orbitals of transition metals, then the projection is automatically 
onto this manifold. Namely, already at the COHSEX level, we expect to recover the benefits 
of an LDA+U type treatment. This is no surprise since a static and localized approximation 
to the GW self-energy yields the kernel of the LDA+U approach [TJ. As an added benefit, 
a COHSEX type approach should automatically include a properly screened U parame- 
ter, as opposed to U being an externally chosen parameter in the usual LDA+U method. 
The dynamic formulae of Eqs. (44) and (45) extend these results to include the frequency 
dependence of the screening. What is missing from these approximate results are the con- 
tributions to the RPA correlation energy which are physically distinct from those stemming 
from screened interband transitions: these are the neglected residue contributions in Eq. (43) 
from the poles from the screening modes themselves. Unfortunately, at present we are un- 
able to provide simple physical examples of situations where these neglected contributions 
play the dominant role and our key approximation to fail for basic physical reasons. This 
question is a subject of present investigation. 



11. NUMERICAL TESTS: ATOMS 

In this section, we describe numerical results on atomic systems. The aim of this section 
is not to present an exhaustive and comprehensive treatment. That is the subject of a future 
investigation. Rather, the main aim is to demonstrate the numerical efficacy of the plasmon 



formula Eq. (36) for the GPU- RPA correlation energy and to test the main approximation 



of high frequency screening used in the approximate results of Eqs. (44) and (45). 

Our atomic code uses a standard non-relativistic approach. The atomic eigenfunctions 
are assumed to take a spherical form: for an eigenstate with spin index a — ±1, the spatial 
part is R n i a (r)Yi m (9,(f)). The radial part R n i(r) is represented numerically on a radial grid 
of exponentially-spaced points El] • The spherical harmonics Y\ m are included up to a 
some maximum angular momentum l max , typically l max = 3 below. When building up the 
one-particle density matrix p a for spin a, 



n,l,m 



Im \ 



we allow for m-dependent state fillings f n ima £ {0, 1}. The local spin-density approximation 
(LSDA) |62j or unrestricted Hartree-Fock (HF) [65] equations are solved by minimizing the 
appropriate total energy over arbitrary occupied orthonormal radial functions. We have 
tested our code with available high quality LSDA [64J and HF [66] data to ensure agreement 
to at least one part in 10 8 in total energies. When computing various correlation energies 
below, we make a spherical approximation which is consistent with assuming spherical eigen- 
functions. This means that when we compute the polarizability (P or P t ) and the screened 
interaction (W or W t ), we assume the fillings have no m-dependence (i.e., shell occupancies 
fnicr instead of orbital occupancies f n i mcT ): different m-channels will not interact so that 
taking both I and m as good quantum numbers is a self-consistent assumption. However, 
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when using matrix elements of W or Wt to compute a contribution to the energy such as in 
Eq. (44), we do include the full m-dependence of the occupancies f n i ma - 



We begin by showing direct numerical evidence for the unboundedness of the GW-RPA 
correlation energy $Jf w that was proven in Section [8] We consider the case of the boron 
atom with configuration ls 2 2s 2 2p 1 . The occupied single-particle states and energies are 
found self-consistently within LSDA for a radial grid extending to r max = 10 Bohr radii. 
After finding the self-consistent potential for the LSDA ground state, we generate the lowest 
300 states for each angular momentum and spin channel. We then compute and tabulate 



the Coulomb matrix elements in the RPA/Casida Eq. (36) within this basis of transitions 



and find the correlation energy of Eq. (38) through direct diagonalization. Starting with this 



information, we then scale all transition energies A t uniformly by a factor A where < A < 1 



so that A t — > AA t . We use the scaled transition energies in the RPA/Casida Eq. (36) to find 



the corresponding plasma frequencies and then use Eq. (38) to compute as a function 
of A. Figure [2] shows the correlation energy as a function of A. As is evident, the correlation 
energy becomes very negative and unphysically large in magnitude. Furthermore, the curve 
is monotonic in A: starting with the reasonable A = 1 energy, one can drive the correlation 
energy to arbitrarily negative values along a continuous path with no extremum along this 
coordinate. As discussed in Section [9l the approximate forms also have this unphysical 



behavior, and the Figure shows the example of the static approximation of Eq. (46). 



We now move onto more pragmatic questions, the first being the convergence properties 



of the plasmon form of Eq. (38). As discussed above, the plasmon form for the GW-RPA 



correlation energy is more physically transparent than the standard integral expression of 



Eq. (23). In addition, it tuns out that it converges much more rapidly. This is demonstrated 
for the cases of the helium atom (Is 2 configuration) and boron atom (ls 2 2s 2 2p 1 configura- 
tion) in Table [TJ To generate the data in this table, an atomic LSDA calculation is run to 
self-consistency for each atom for a radial grid of r max = 7 Bohr radii for He and r max = 10 
for B, and the lowest 300 eigenstates and eigenenergies of the Kohn-Sham Hamiltonian are 
computed and tabulated. The N max lowest-energy eigenstates are then used to evaluate the 
correlation energies, and convergence is monitored by increasing N max . 



For the plasmon form, the RPA/Casida Eq. (36) is solved within the basis of transitions 
generated by the N max single-particle states via direct diagonalization. For the integral form, 
the computations are more involved. First, for numerical stability, the integral along the 
real axis is changed to along the imaginary uj axis (renamed /3) via a Wick rotation. Using 
the identity tr In A = In det A for a matrix A, the integral is then 

®c W = 4£<W> + \ /"f hidet^/3). (53) 



t 



— oo 



The matrix determinant is computed within the subspace spanned by the single-particle 
states. Second, the upper limit is reduced to a finite value (3 max and the integral is discretized 
with spacing A/3. As we desire the integral in the limits A/3 — > and (3 max — > oo, two 
extrapolations are performed. For fixed f3 max , the integral is evaluated for a series of A/3, 
and Richardson extrapolation is performed to A/3 = 0. Next, these extrapolated values 
are themselves Richardson extrapolated to (3 max = oo by noting that for large (3, e(i(3) = 
I + A(3~ 2 + 0(/3~ 4 ) so that the neglected integral from (3 max to oo is proportional to /3~„ x 
to leading order. These extrapolation procedures are straightforward and unproblematic 



29 




Boron Atom Correlation Energy 



□ d □ □ ^ 



-0--0--0 



GO Static 
QHRPA 



0.4 0.6 
X 



0.8 



FIG. 2: Correlation energy <3? c as a function of the scaling A of the transition energies for the boron 
atom. Single-particle energies are from a ground-state LSDA calculation. The blue squares with 



solid line are the exact RPA correlation energies from Eq. ( 38 ) when using scaled transition energies 



in the RPA/Casida Eq. (36). The green squares with dashed line are correlation energies from the 



static approximation of Eq. (46) using the scaled transition energies for computing the screening. 



The dashed or solid lines are guides for the eye. Clearly, the energy decreases monotonically with 
decreasing A to large and unphysical negative values with no extrema along this path. 



because the integrand is smooth as a function of /3. The doubly extrapolated values are 
listed in Table [H 

Clearly, the plasmon form has superior convergence properties when compared to the 
integral form, and, additionally, does not require any particular extrapolation or discretiza- 
tion procedures. Helium presents a rather simple case: more complex atoms (or molecules) 
will require progressively larger basis sets to converge the integrals and the plasmon form 
should prove even more useful in practice. A hint of this is provided by comparing He to B 
in Table [Tj For B, the integral form of the correlation energy converges more slowly, and the 
calculations at N max = 300 were already quite demanding in terms of time (and patience). 
The plasmon form, again, shows rapid convergence versus N max . 

We now examine the accuracy of the approximate forms derived in Sections [9] and 10 
above. In Tables [IT] and HI we present these correlation energies for the helium and boron 
atoms using LSDA or HF wave functions and eigenenergies. These are to be compared to 
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TABLE I: Comparison of convergence of GW-RPA correlation energies for the He and B atoms. 



The RPA correlation energy is computed using the standard integral of Eq. ( 53 ) (third column) or 



the plasmon form of Eq. (38) (fourth column). All energies are in Hartrees. N max is the number 
of single-particle state included in the calculations for each angular momentum and spin channel. 
The single particles energies and wave functions are from the LSD A. 



Atom 




integral 


plasmon 


He 


25 


-0.1038 


-0.0804 


He 


50 


-0.0916 


-0.0805 


He 


100 


-0.0849 


-0.0806 


He 


150 


-0.0830 


-0.0806 


He 


200 


-0.0819 


-0.0806 


B 


50 


-0.4126 


-0.2129 


B 


100 


-0.3600 


-0.2171 


B 


200 


-0.2844 


-0.2175 


B 


300 


-0.2584 


-0.2175 



the exact RPA energy of Eq. (38) in the second to last row in each table. For completeness, 
we also include the other energy terms to show their relative importance and their variation 
with single-particle theory, discussed further below. The Hartree energy reported in the table 
is the one based on the actual, non-spherical electron density {i.e. based on m-dependent 
occupancies fnima as opposed to the spherical density in most DFT atomic calculations). 
We begin by considering our main approximate form of Eq. (45). As we can see from 



comparing to the exact RPA energy, and especially when comparing to the static versions, 



the basic approximation underlying Eqs. (44) and (45) is a relatively good one: the absolute 
correlation energy of Eq. (45) differs by at most 0.1 Ha from the exact RPA one. This shows 
that the fundamental approximation of assuming high-frequency screening is reasonable even 
in atoms. For solids and extended systems where the Coulomb interaction shows true long- 
ranged behavior (as opposed to atoms), plasma modes are of higher energies than interband 
energies and the situation should be further improved. 



As explained above, the approximate form of Eq. (45) is as computationally expensive to 
calculate as the exact RPA plasmon form but is additionally very difficult to converge. To 
obtain the values in the tables, we had to perform the following steps simultaneously: (i) 
increase the size of the radial axis to make for a denser continuum, (ii) increase the number 
of one-particle states entering the calculation to ensure a fixed level of convergence with 
increasing radial axis, and (hi) exclude contributions to W t (A t ) from transitions t' that were 
within a small energy window 5 of A t (i.e., \A t > — A t | < 5) while sending 5 — > 0. 

Moving on to the approximations that assume static screening, we see that they overes- 
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TABLE II: Various components of the total energy of the helium atom. Results are reported for 
single-particle wave functions and eigenenergies coming from self-consistent LSDA or Hartree-Fock 
(HF) calculations. Energies are in Hartree. Except for the (*) values, all calculations use a radial 
grid of size r max = 18 Bohr radii and N max = 257 single-particle states for each angular momentum 
channel. The starred (*) values are difficult to converge, and the values in the table are uncertain 
to within ±0.001: see text for details. For reference, the exact (CI) non-relativistic correlation 
energy is also shown as the last entry. 



Energy component 



Equation LSDA 



HF Difference 



(22) 



Kinetic 
Electron-ion 
Hartree 
Fock exchange 

$ c : W(0) instead of W t (0) (|50|) 
<]? c : square root approx. 
<]? c : static approx. 
<E> C : main approx. 
$ c : RPA 

$ c : CI (exact) Ref. [67] 



(49) 



(46) 



(45) 



(38) 



2.768 
-6.626 

1.996 
-0.998 
-0.318 
-0.311 
-0.313 
-0.060* 
-0.081 

-0.0420 



2.862 
-6.749 
2.052 
-1.026 
-0.255 
-0.248 
-0.250 
0.048* 
-0.064 



1.0% 
1.0% 
1.0% 
1.0% 
25% 
25% 
25% 
25% 
27% 



timate the importance of correlations and give too negative values uniformly. We discuss 
the reason for this in the next paragraph, but in the meantime we see that once the static 
approximation is made, the various forms for the correlation energy are quite similar. For 



example, the difference between the square root form of Eq. (46) and its series expansion in 



Eq. (49) is small. The smallness of the differences simply means that the ratio (W(0))/ A is 



small: e.g., for boron we find the largest value of the ratio (W(0)}/A is 0.2 and is achieved 
for the 2s-2p transition. Given how far all these static correlation energies are from the 



dynamic answer of Eq. (45) or the exact RPA answer of Eq. (38), one can take all these 



static approximations to be basically of equal accuracy. 

The reason the static approximations overestimate the magnitude of the correlation en- 
ergy is easy to understand. Figure [3] shows how three representative correlation energies 



converge for the case of atomic boron: the static formula of Eq. (46), the dynamic formula 



of Eq. (45), and the exact RPA formula of Eq. (38). What is shown is the correlation energy 



contributions summed up to a given transition energy (or plasmon energy for the exact case). 
In all cases, we see that the contributions to the correlations are small, then become large at 
a certain energy, and then become smaller again. In atomic boron at LSDA level, there are 
two physically important transitions: the dominant 2s-2p at 0.21 Ha and then the weaker 
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TABLE III: Various components of the total energy of the boron atom. Results are reported for 
single-particle eigenfunctions and eigenenergies coming from self-consistent LSDA or HF calcula- 
tions. All energies are in Hartree. Except for the (*) values, all calculations use a radial grid of size 
r max = 40 Bohr radii and N max = 400 single-particle states for each angular momentum channel. 
The starred (*) values are difficult to converge, and the values in the table are actually uncertain to 
within ±0.01: see text for details. For reference, the exact (CI) non-relativistic correlation energy 
is also shown as the last entry. 



Energy component 



Equation LSDA 



HF Difference 



(22) 



Kinetic 
Electron-ion 
Hartree 
Fock exchange 

$ c : W(0) instead of W t (0) (|50|) 
<]? c : square root approx. 
<]? c : static approx. 
<E» C : main approx. 
$ c : RPA 

$ c : CI (exact) Ref. 



(49) 



(46) 



(45) 



(38) 



24.173 
-56.520 
11.534 
-3.712 
-0.869 
-0.820 
-0.835 
-0.30* 
-0.217 

-0.125 



24.530 
-56.900 
11.590 
-3.749 
-0.717 
-0.680 
-0.692 
-0.11* 
-0.171 



1.5% 
0.67% 
0.49% 
1.0% 
21% 
21% 
20% 
270% 
27% 



ls-2p at 6.4 Ha. As the transition energy sweeps through each atomic transition, we see 
large contributions to the correlation energies at first, but then the transition is "exhausted" 
and the contributions become small again. 

The obvious difference between the three curves is that for the static case all the con- 
tributions are always negative, tend to be large in magnitude, and keep adding up to yield 
a large negative value. On the other hand, for the dynamic approximation and the exact 
result, the initial low-energy contributions are large and negative, may then become slightly 
positive, but rapidly become small in magnitude. To understand this difference, we write 
the modified screened interaction as W t (u) = ^(u;) -1 ^ in terms of a modified dielectric 



function e t (u). Solving the RPA/Casida Eq. (36) with transition t missing gives us a set 
modified plasma modes with energies u p and mode functions \p). Based on Eq. (35), we 
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FIG. 3: Convergence of correlation energies for the boron atom. Both plots show the cumulative 
sum of the contributions to the correlation energies (vertical) for transitions and/or plasma modes 
up to some given energy At (horizontal). LSDA wave functions and eigenenergies are used with 
a radial grid of size r max = 40 Bohr radii and N max = 400 eigenstates. The lowest solid black 



curve is for the static approximation of Eq. ( 46 ) , the middle red dashed curve is for the dynamic 



approximation of Eq. (45), and the uppermost green dotted curve is for the exact RPA formula 
of Eq. (38). The plots show the same data with the only difference being a linear (top plot) or a 
logarithmic scale (bottom plot) of the horizontal axis. The dominant transitions in LSDA boron 
are the 2s-2p at 0.21 Ha and ls-2p at 6.4 Ha, both visible in the lower plot as energies where the 
correlation contributions have a sudden jump. 

have the dynamic and static matrix elements 

p 1 p 

2u p \(t\V\p)\* 



(t\w t (o)\t) = (t\v\t)-Yl 



P v 

Clearly, the static formula always gives negative correlation contributions that become small 
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only when the matrix elements (i|V|p) become small. On the other hand, the dynamic case 
has a denominator that changes sign from positive to negative as A t sweeps through the 
plasma energies, and the denominator itself gets large when A f gets large. The dynamic 
screening behavior of the plasmons is missing in the static formula which assumes that 
no matter what the transition energy, the plasmons can screen it adiabatically. This is 
obviously erroneous for transitions where A t > Q p . Therefore, the static formula can be 
improved by adding some dynamic behavior in the screening even if done approximately. 
This is an example of what was meant in Section 10 regarding the use of plasmon-pole 
models or model dielectric functions to improve the static COHSEX. 

We end this section with some observations on the results in Tables [Til and II I II and their 
dependence on the single-particle theory. The parts of the total energy that depend only 
on the single-particle orbitals, i.e. the kinetic, electron-ion, Hartree, and Fock exchange 
energies, depend weakly on the choice of LSDA versus HF single-particle orbitals, changing 
at most ~1%. This is not surprising a posteriori as visual comparison of the radial functions 
show small differences. However, the correlation energies depend more strongly on the 
choice of single-particle theory, and this is due to the relatively large differences between 
the single-particle energies. The HF-based correlation energies are smaller than the LSDA 
simply because HF transition energies are larger than LSDA: e.g., for atomic boron the 
important 2s-2p transition is at 0.210 Ha in LSDA but at 0.232 Ha in HF; larger transition 
energies mean weaker screening and thus weaker correlation. As per Section |8j it matters 
greatly whether the Green's function Go is generated by a local or non-local potential. 

The true many-body Green's function G obeys Dyson's Eq. Q and is thus generated by a 
non-local (and dynamic) self-energy S xc (x, x', u). A priori, we would expect that a static but 
non-local potential Uq(x,x') should generate a decent non-interacting Go which should be 
close to the true G, certainly closer than one generated by a static local potential. However, 
as discussed in Section [8j choosing the "best" Uo via an optimization of the total energy 
functional is problematic unless constraints are imposed on Uq. The outstanding theoretical 
problem is what constraints to impose and which ones are the "best" ones. This obviously 
involves creating and justifying metrics that tells us how good each choice of constraints will 
be in practice. 



12. SUMMARY AND OUTLOOK 

Our work has focused on the correlation energy functional with Luttinger-Ward theory 
(specifically the Klein functional) within the GW-KPA approximation. The Green's func- 
tions used in the energy functional are of non-interacting form. The main findings in this 
work are threefold. First, we present the exact rewriting of the Giy-RPA correlation en- 



ergy functional in Eq. (38) in terms of differences between plasma and interband energies. 
This form is directly amenable to computation, shows good convergence properties in the 
atomic tests, and has prospects for having its computational scaling improved by use of ma- 
trix square root algorithms. Second, we describe the approximate rewriting of the GH^-RPA 



correlation energy functional in Eqs. (44) and (45) where the correlation energy is written as 



a sum of screened interband transition contributions; the main approximation is to assume 
that the dominant screening dynamics are much faster than the key interband dynamics; 
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atomic tests show that the approximation is good numerically. These approximate forms 
then lead to a ladder of approximations where the COHSEX is the simplest one possible. 
Third, we show, analytically and with numerical examples for atoms, that if one restricts the 
Green's function to be of non-interacting form generated by Hermitian non-local potentials, 
the GW-RPA correlation energy has no lower bound over this set of Green's functions; nor 
does it have an extremum. 

Going forward, the exact rewriting and its good convergence properties — coupled with 
algorithmic development — should pave the wave for wider application of the GW-RPA 
correlation functional to materials systems. The multitude of approximate forms we present 
here will hopefully broaden and improve the types of approximate self-energies used in 
self-consistent band structure methods that go beyond the usual LDA or GGA treatments. 
Simultaneously, having explicit energy functionals on hand means one can construct the self- 
energies in a variational manner: the solution of the self-consistent equation for a particular 
self-energy will optimize the total energy functional from which it was derived. However, 
the result on the unboundedness of the GW-RPA correlation energy over the space of all 
non-interacting Green's functions means that applying these Green's function approaches is 
not yet a rote exercise in applying standard unconstrained optimization algorithms. Rather, 
the choice of non-local potential that generates the non-interacting Green's function must 
be constrained in some manner so as to avoid the pathological negative infinite correlation 
energy and to produce extrema. The theoretical question is then to understand which 
constraints succeed and also to devise metrics to compare their quality and accuracy. 
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